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d We calculate cumulants of fluctuations of net-baryon number, net-electric charge and net- 

strangeness, in the framework of lattice regularized QCD. We use a highly improved staggered 

00 

quark (HISQ) action on lattices with temporal extent of A^^ =6,8 and 12 and almost physical 
quark masses. By means of a Taylor expansion in various chemical potentials and under de- 
manding both strangeness neutrality as well as the correct isospin asymmetry, we evaluate these 
fluctuations at conditions met in heavy ion collisions. Cumulants of net-electric charge fluctu- 
ations can, in principle, also be measured in heavy ion experiments. We therefore propose a 
method to extract freeze-out parameters, such as the freeze-out temperature and baryon chemical 
potential, based on a comparison of lattice results and experimental measurements of two differ- 
ent ratios of net-electric charge cumulants. As this method involves only (lattice) QCD results and 

>- 

\^ experimental measurements it is model independent and, for the purpose of extracting freeze-out 

04 parameters, does not require any input from hadron resonance gas model calculations. 
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1. Introduction 

To obtain a detailed and quantitative understanding of tiie QCD phase diagram is one of the 
most important and outstanding problems in high energy physics. The analysis of the phase dia- 
gram is naturally connected with an investigation of the properties of strongly interacting matter at 
high temperature and densities and the nature of the QCD phase transition. One expects that the 
nature of the transition depends on the quark masses and chemical potentials. A generic phase dia- 
gram based on model calculations and model independet arguments is shown in Fig. 1 [1]. One of 
the most interesting question is to clarify whether there exists a critical point in the phase diagram 
at physical quark masses or not. 



Figure 1: Generic phase diagram of QCD, based on model calculations and model independent symmetry 
arguments. Also indicated are the regions in the phase diagram where we are able to obtain results on 
fluctuation observables from lattice QCD and experiments, respectively. 

Recently large efforts are made in probing the QCD phase diagram with heavy ion collision. 
The RHIC beam energy scan program (BBS) aims at a systematic scan of the QCD phase diagram. 
Through a variation of the center of mass energy of the two colliding ions, the fireball made of 
quark gluon plasma (QGP) is generated under different initial conditions, as e.g., different initial 
temperatures and net-baryon densities. After its formation, it evolves on isentropic (hydrodynamic) 
trajectories in the phase diagram, i.e. it expands and cools. It eventually goes through the QCD 
transition and hadronizes. The point from whereon the abundance of the hadronic particle species 
do no longer change is called chemical freeze-out. The particle yields realized here are then fi- 
nally measured by various detectors, only modified by particle decays. One expects that under the 
variation of the collision energy {y/sAA), the freeze-out happens along a one-dimensional curve in 
the phase diagram that can be parametrized by yJs^A [2]- The freeze-out curve is also shown in 
Fig. 1. One way to calculate the freeze-out parameters such as temperature, chemical potentials 
and volume is by performing a least-square-fit of the Hadron Resonance Gas (HRG) model to the 
measured particle yields. This procedure has been very successful in the past [3]. 
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Theoretically best understood is the QCD phase diagram at vanishing net-baryon density, or 
equivalently zero baryon chemical potential. Here lattice regularized calculations of QCD are 
feasible and have reached a precision level that enables controlled continuum extrapolated results 
with physical quark masses. Also this region is indicated in Fig. 1. Unfortunately, lattice QCD 
calculations are not possible at non-vanishing net-baryon densities by means of standard Monte 
Carlo methods, due to the notorious sign problem. Nevertheless, there are possibilities to extract 
lattice QCD results at small but nonzero baryon chemical potentials. A straightforward approach is, 
e.g., given by expanding the observables in a Taylor series around zero baryon chemical potential 
(/Ifi = 0) [4]. A method to calculate the freeze-out curve by a direct comparison of lattice results 
with experimental measurements has been recently proposed by us [5, 6, 7, 8]. It is based on 
a Taylor expansion of ratios of cumulants of net-charge fluctuations. Since it involves only first 
principle lattice QCD calculations it is (HRG) model independent. 



2. Cumulants of conserved charges 

We start with considering an expansion of the QCD partition functions - or rather the pressure 
in dimensionless units - in terms of the net-baryon, net-electric charge and net-strangeness chem- 
ical potentials, also in dimensionless units (/ts = ^B/T,(iQ = ^q/T^fis = l^s/T) which is given 
by 

^ = ^ i'^^ = £ ^xfZ Ai al ■ (2. 1) 

We calculate the temperature dependent coefficients xfj^^Q the lattice at vanishing chemical po- 
tentials. On the lattice we obtain them as derivatives of the partition function (generalized suscep- 
tibilities) exactly as given by the definition of the Taylor series 



. . (2.2) 

jl=o 



Interestingly, these coefficients also define cumulants of the corresponding net-charge fluctuations 
{Nb,Nq,Ns), for the diagonal coefficients one finds 

iVT')-x! = {(SNxf), (2.3) 
{VT')-x! = {i6Nx)')-3{{8Nx)'y, (2.4) 
{VT')-xi = {idNxf)-l5(^{6Nx)'){idNxf)+30(^{dNx)')\ (2.5) 

with X = B,Q,S and 5Nx = Nx — {Nx). Similar relations can also be found for the off-diagonal 
coefficients. The cumulants of charge fluctuations can in principal be measured in heavy ion ex- 
periments. However, as experiments have difficulties to trigger on neutrons and various strange 
particles, baryon number and strangeness fluctuations seem to be unfeasible. On the other hand, 
experiments can measure proton number fluctuations [9], which can not be calculated in QCD since 
the proton number is not a conserved quantity. Unless one finds a way to relate proton number to 
baryon number fluctuations [10, 11] the only remaining set of cumulants that can be obtained in 



3 



Conserved charge fluctuations with HISQ fermions 



Christian Schmidt 



both experiments and lattice QCD calculations are the net-electric charge fluctuations. Neverthe- 
less, we have to keep in mind, that the measurement of the fluctuations of a conserved charge is a 
difficult experimental task, as it involves a delicate tuning of the rapidity window that is accepted in 
the detectors [12]. Choosing the acceptance range too large, will result in a suppressed signal as the 
global charge conservations forces the fluctuations to vanish in the limit of a total {An) acceptance, 
on the other hand, an acceptance that is too small will loose the sensitivity to the relevant physics 
and/or not survive hadronization and the hadronic phase. 

In lattice QCD calculations electric charge fluctuations have systematic errors that need to be 
controlled. Although statistical errors are much smaller than in the baryon number fluctuations, 
one has to work on very fine lattices in order to control the continuum limit systematically. This is 
due to the fact that large contributions to the electric charge fluctuations stem from the light pion 
sector. The pion spectrum is distorted on the lattice, when one uses the popular staggered fermion 
formulation of the QCD action. For reasons of the numerical costs, staggered fermions are the 
most common type of fermions for thermodynamic calculations on the lattice. In the following we 
will use the highly improved staggered fermion (HISQ) action [14], with reduces the deformation 
of the pion sector as much as possible. 

In Fig. 2 we show our preliminary results of the cumulants (expansion coefficients), obtained 
with the HISQ action. The calculations have been done with 2 light and one heavier quark flavor 
and almost physical quark masses, i.e. the strange quark mass has been fixed to its physical value, 
whereas the light quark mass was taken to be 1/20 of the strange quark mass. In general we find 
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Figure 2: Second and fourth order cumulants of net-baryon number fluctuations (left column), net-electric 
charge fluctuations (middle column) and net-strangeness fluctuations (right column). Different symbols 
denote different lattice spacings: A^^ = 6 (triangles) and 8 (circles). The vertical bar indicates the QCD 
transition temperature as obtained in [15]. 

little dependence on the lattice cutoff, which is here controlled by the number of lattice points in 
(Euclidean) temporal direction (A^t)- In case of the 2"'' order fluctuations the HotQCD Collabo- 
ration has performed a continuum extrapolation based on A^,- = 6,8 and 12 lattices [16]. Good 
agreement with HRG model results has been found for temperatures below T < 150 MeV. 

The structure of the cumulants is found to be consistent with the expected QCD critical be- 
havior connected to the 0(4) critical point in the massless limit of the two light quark flavors (see 
Fig. 1). From an analysis of the free energy in terms of universal scaling fields one would expect 
the singular part of the cumulants to behave as ~ |f p^"^" [13], at least for X = B and Q, where 
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t is the reduced temperature and a the critical exponent of the specific heat. As a is small and 
negative (a —0.15) the 4"^ order fluctuations should develop a cusp in the chiral limit whereas 
the 6* order fluctuations are the first to diverge. 

3. Electric charge and strangeness chemical potentials 

In order to resemble the conditions met in heavy ion collisions as closely as possible, we 
demand strangeness neutrality {{Ns) = 0) and the correct isospin asymmetry (r = {Nq) / {Nb))- 
These two conditions can be realized by choosing the free parameter {Xq and jXs accordingly. By 
expanding both the conditions as well as the chemical potentials in terms of {Xb, 



fls = si{T)flB + S3{T)fll + Ae = ^1 {T)flB + c]3{T)fL^ + ml 



r.3 



a5\ 



(3.1) 



we can solve for the coefficients in Eq. 3. 1 order by order and thus enforce the conditions up 
to arbitrary order in /tfi. Or results for the leading order (LO) and next to leading order (NLO) 
coefficients of the /tg and jXs series are shown in Fig. 3 (left and middle). The upper panel shows the 
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Figure 3: The leading and next-to-leading order expansion coefficients of the strangeness (left) and the 
negative of the electric charge chemical potentials (middle) versus temperature for r — 0.4. For si and qi the 
LO-bands show results for the continuum extrapolation. For st, and ^3 we give an estimate for continuum 
results (NLO bands) based on spline interpolations of the Nx = % data. Dashed lines at low temperature are 
from the HRG model and at high temperature from a massless, 3-flavor quark gas. The right hand panel 
shows NLO results for jJ-s/fiB and fig /lis as function of jis for three values of the temperature. 



leading order (LO), whereas the lower panels show the ratio of NLO to LO coefficients. The band 
in the upper panels indicates the continuum extrapolation based on the A^^ = 6, 8 and 12 data, the 
band in the lower panels is a continuum estimate based on spline interpolations of the = S data. 
We find that the NLO contributions are negligible in the high temperature region and below 10% in 
the temperature interval relevant for the analysis of freeze-out conditions, i.e., T ^ (160 it 10) MeV. 
In fact, in this temperature range the leading order lattice QCD results deviate from HRG model 
calculations expanded to the same order by less than 15%. Note that one can also investigate the 
convergence properties of the HRG model itself. In the HRG model the NLO expansion reproduces 
the full HRG result for jig and /tj to better than 1.0% for all values of /ig < 1.3. Altogether, we 
thus expect that the NLO truncated QCD expansion is a good approximation to the complete QCD 
results for fig and fis for ^b < 200 MeV. 
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Our results for the strangeness and electric charge chemical potentials at NLO as function of 
jUg and T are shown in Fig. 3 (right). While pis/^B varies between 0.2 and 0.3 in the interval 150 
MeV < r < 170 MeV, the absolute value of Hq/Hb is an order of magnitude smaller. Both ratios 
are almost constant for ^Ug < 200 MeV, which is consistent with HRG model calculations. 

4. Comparison with the experiment 

We will now construct the observables that we want to compare with the experiment in order 
to determine the remaining freeze-out parameters, which are the freeze-out temperature (T-f), the 
freeze-out baryon chemical potential (/i^) and the freeze-out volume (V^). The latter one can 
be easily eliminated by considering ratios of cumulants as should be apparent from Eqs. 2.2 - 
2.5. We ai^e thus left with two freeze-out parameters (r^,jUg), for which we need two independent 
observables to match with the experiment. As already discussed in Sec. 2, from the set of cumulants 
here considered, only the net-electric charge fluctuations can be determined on the lattice as well 
as measured in experiments. We therefore propose the following two ratios of net-electric chai^ge 
fluctuations for the comparison with the experiment. 



Here we expressed the cumulant ratios also in terms of the mean value (M), the variance (a) and 
the skewness (S), which characterize the shape of the net-electric charge distribution. We have 
again expanded these two quantities in terms of fls and determined the LO and NLO contributions 
in that series (the remaining chemical potentials jJ-g and /I5 have been fixed as described in Sec. 3). 
The two ratios defined in Eq. 4. 1 and 4.2 represent the most simple choice as they only involve 
the evaluation of cumulants up to the 2"'* and 4* order, respectively, at LO and up to the 4"^ and 
6 order, respectively, at NLO. Morover, they are complementary in the sense thatTefj (Eq. 4.1) IS 
an odd function of /Xg, which leads to a distinct sensitivity with respect to the freeze-out chemical 
potential (jU^), whereas R^i (Eq. 4.2) is an even function of jXs, starting with a constant, which 
results in a much more pronounced sensitivity to the freeze-out temperature (T^). 

In Fig. 4 (left) we show our results on the LO and NLO expansion coefficients of 7?^2- The 
bands on the lower and upper panel have the same meaning as in Fig. 3 (left). We find that the 
NLO corrections to Rfj are below 10%, which makes the LO result a good approximation for a 
large range of /tg. Systematic errors arising from the truncation of the Taylor series for i?^2 at 
NLO may again be estimated by comparing the full result in the HRG model calculation with 
the corresponding truncated results. Here we find for T = (160 it 10) MeV and fLs < 1-3 that 
the difference is less than 1.0%. Moreover, we estimated that taste violation effects in the NLO 
calculation lead to systematic eiTors that are at most 5% and thus will be negligible in /?^2 • Taylor 
series truncated at NLO are thus expected to give a good approximation to the full result for a wide 
range of baryon chemical potentials. 

In Fig. 4 (middle) we show the full /ig and T dependence of /?^2' including the NLO con- 
tribution. Obviously the ratio /?^2 shows a strong sensitivity on /Xg but varies little with T in the 




(4.1) 




= R^f>+R^l^fLl + ff{fl',). 



(4.2) 
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Figure 4: The left panel shows LO (top) and NLO (bottom) expansion coefficients of for r — 0.4. The 
bands and lines are as in Fig. 3(left).On the middle panel we plot versus P-b/T, including the NLO 
contribution, for three values of the temperature. R^^ versus temperature is shown on the right panel for 
/Ib = 0. The wider band on the data set for Ni = % shows an estimate of the magnitude of NLO corrections. 



temperature range T = (160 it 10) MeV. For the determination of {T^ a second, complimen- 
tary information is needed. To this end we use the ratio R^^, which is strongly dependent on T but 
receives corrections only at ^(/t|). The leading order result for this ratio is shown in Fig. 4 (right). 
Apparently this ratio shows a characteristic temperature dependence for T > 155 MeV that is quite 
different from that of HRG model calculations. The NLO correction to this ratio vanishes in the 
high temperature limit and at low T the HRG model also suggests small corrections. In fact, in 
the HRG model the LO contributions to 7?^j differ by less than 2% from the exact results on the 
freeze-out curve for pLs < 200 MeV. The broader band in Fig. 4 (right) indicates an estimate of the 
NLO contribution at /ts = 1 from our Nt = S calculations. 

We now are in the position to extract jU^ and Tf from /?^2 ^3^1 which eventually will be 
measured in the beam energy scan at RHIC [17, 18]. A large value for R^^, i.e. R^^ ~ 2 would 
suggest a low freeze-out temperature T < 155 MeV, while a value R^^ ~ 1 would suggest a large 
freeze-out temperature, T ~ 170 MeV. A value of R^^ ~ 1.5 would correspond to T ~ 160 MeV. 
A measurement of R^^ thus suffices to determine the freeze-out temperature. In the HRG model 
parametrization of the freeze-out curve [2] the favorite value for T-^ in the beam energy range 
200 GeV > y/sXX > 39 GeV indeed varies by less than 2 MeV and is about 165 MeV. At this 
temperature the values for R^^ calculated in the HRG model and in QCD differ quite a bit, as 
is obvious from Fig. 4 (right). While R^y ~ 2 in the HRG model, one finds ~ 1.2 in QCD 
at r = 165 MeV. Values close to the HRG value are compatible with QCD calculations only for 
r < 157 MeV. We thus expect to either find freeze-out temperatures that are about 5% below HRG 
model results or values for Rf^ that are significantly smaller than the HRG value. A measurement 
of this cumulant ratio at RHIC thus will allow to determine and probe the consistency with 
HRG model predictions. 

For any of these temperature values a comparison of an experimental value for R^2 with Fig.4 
(middle) will allow to determine /Xg. To be specific, at T = 160 MeV we expect to find jx^ = 
(20 - 30) MeV, if /?f2 lies in the range 0.012 - 0.020, = (50 - 70) MeV for 0.032 < R^^ ^ 0-045 
and njj = (80 — 120) MeV for 0.05 < /^fj < 0.08. These parameter ranges are expected [2] to cover 
the regions relevant for RHIC beam energies ^/sAA = 200 GeV, 62.4 GeV and 39 GeV, respectively. 
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As is evident from Fig. 4 (middle) the values for yi^ will shift to smaller (larger) values when 
turns out to be larger (smaller) than 160 MeV. A more refined analysis of (r^,/x^) will become 
possible, once the ratios ^^id R^^ have been measured experimentally. 
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